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Abstract — We present a rigorous and computationally ef- 
ficient method to do a parameter-free analysis of molecular 
wires connected to contacts. The self-consistent field ap- 
proach is coupled with Non-equilibrium Green's Function 
(NEGF) formalism to describe electronic transport under 
an applied bias. Standard quantum chemistry software is 
used to calculate the self-consistent field using density func- 
tional theory (DFT). Such close coupling to standard quan- 
tum chemistry software not only makes the procedure sim- 
ple to implement but also makes the relation between the 
I-V characteristics and the chemistry of the molecule more 
obvious. We use our method to interpolate between two ex- 
treme examples of transport through a molecular wire con- 
nected to gold (111) contacts: band conduction in a metallic 
(gold) nanowire, and resonant conduction through broad- 
ened, quasidiscrete levels of a phenyl dithiol molecule. We 
obtain several quantities of interest like I-V characteristic, 
electron density and voltage drop along the molecule. 

Keywords — Molecular Electronics, First-Principles, DFT, 
NEGF 

PACS- 85.65.+h, 73.23.-b, 31.15.Ar 

I. Introduction 

There is much current interest in molecular electron- 
\ ics due to the recent success in measuring the I-V char- 
acteristics of individual or small groups of molecules Q, 
H> Hi , @> @, @, 0: §■ The measured resistances ex- 
hibit a wide range of values. For example, n-alkane chains 
(CH3 — (C-ff 2 )n-i) have large gaps (6 eV or greater) be- 
tween the highest occupied molecular orbital (HOMO) and 
the lowest unoccupied molecular orbital (LUMO) and act 
like strong insulators ||, |1(J while gold nanowires (or 
quantum point contacts) have zero gap (or a continuous 
density of states near the Fermi energy) and exhibit novel 
one-dimensional metallic conduction characteristics JlTJ ] , 
p2[ . Similar one-dimensional metallic conduction (but on 
a much larger length scale of the order of micrometers) has 
been observed in carbon nanotube s |13[ ]. Then there are 
biological molecules like the DNA |L4|], fig] whose electri- 
cal conduction characteristics are currently the subject of 
much debate |ll| . Interesting device functionalities such as 
transistor action have also been recently reported ||l7| . 

Understanding the correlation between the chemical 
and electronic properties of a wide class of molecules is 
a first step towards developing a 'bottom-up' molecule- 
based technology. Semi-empirical theories [Q, |Q, po| , 
@, (22), ||, @, ||, gf| have been used to qualita- 
tively study electronic transport in molecules. Some first- 
principles methods have also been developed |^7j , Q , Q , 
[j30| , |3ll] . Typically the first-principles methods are either 
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Fig. 1 

Schematic of a molecule (phenyl dithiol) connected to two 
semi-infinite gold contacts. the self-energy matrices si and 

£2 exactly ACCOUNT FOR THE CONTACTS AND ARE OF THE SAME SIZE 

AS THE MOLECULAR FOCK MATRIX F. THE SELF-CONSISTENCY 
SCHEME IS SHOWN FOR (a) ISOLATED MOLECULE IN EQUILIBRIUM AND 
(b) CONTACT-MOLECULE- CONTACT SYSTEM UNDER BIAS. 



computationally very expensive or ignore charging effects 
by employing non self-consistent calculations . Further- 
more, the correct geometry and atomicity of the contacts 
is often ignored by employing a jcllium-like model although 
surface effects like bonding and chemisorption are clearly 
very important. The purpose of this paper is to describe 
a straightforward (computationally inexpensive) yet rigor- 
ous and self-consistent procedure for calculating transport 
characteristics while taking into account the effects men- 
tioned above. In contrast to an ab-initio theory simulating 
the large (in principle infinite) open system, we partition 
the problem into a 'device' and a 'contact' subspace such 
that standard quantum chemistry techniques can be em- 
ployed to analyze the electronic structure of a finite-sized 
device subspace, while incorporating the effects of the out- 
side world ('contacts') through appropriately defined self- 
energy matrices and fields. 

The electronic structure of an isolated device is obtained 
in standard quantum chemistry through a self-consistent 
procedure pl| , |3j| shown schematically in Fig [j] a. The 
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process consists of two steps that have to be iterated to 
obtain a self-consistent solution. Step 1: Calculate the 
Fock matrix F for a given density matrix p using a specific 
scheme such as Hartree-Fock (HF) or density functional 
theory (DFT) to obtain the self-consistent field , and Step 
2: Calculate the density matrix from a given self-consistent 
Fock matrix F and total number of electrons N, based on 
the laws of equilibrium statistical mechanics. 

The problem of calculating the current-voltage (I-V) 
characteristics is different from the above self-consistent 
procedure in the following ways : (i) We are dealing with 
an open system having a continuous density of states and 
variable (in principle fractional) number of electrons, rather 
than an isolated molecule with discrete levels and integer 
number of electrons; (ii) The molecule does not necessarily 
remain at equilibrium or even close to equilibrium - two 
volts applied across a short molecule is enough to drive it 
far from equilibrium; (iii) Surface effects like chemisorp- 
tion and bonding with the contacts are expected to play 
a non-trivial role in transport. To calculate the density 
matrix (step 2 above) we therefore need a method based 
on non-equilibrium statistical mechanics applicable to such 
an open system with a continuous density of states. The 
Non-equilibrium Green's Function (NEGF) formalism p4| , 
[ [35t provides us with such a method and that is what we 
use in this paper for step 2. Step 1, however, remains un- 
changed and we use exactly the same procedure as in stan- 
dard quantum chemistry software. The overall procedure 
is shown schematically in Fig [l] b. 

We rigorously partition the contact-molecule-contact 
system (see Fig Q) into a contact subspace and a device 
subspace. This makes the molecular chemistry conceptu- 
ally transparent as well as computationally tractable. The 
contact subspace is treated via a one-time calculation of 
the surface Green's function of the contacts including their 
atomicity and crystalline symmetry. Different molecules 
coupled to the same contacts have different couplings but 
the contact surface Green's function is independent of the 
molecule. Given the surface Green's function and the 
contact-molecule coupling we can describe the chemisorp- 
tion and bonding of the molecule with the infinite contacts 
through self-energy matrices of finite size (equal to that of 
the molecular subspace). The NEGF formalism has clear 
prescriptions (Section [ij) to calculate the non-equilibrium 
density matrix from a knowledge of F, and the elec- 
trochemical potentials in the two contacts, fi% and p2- All 
quantities of interest (electron density, current etc) are then 
calculated from the self-consistently converged density ma- 
trix. 

From a computational viewpoint, the most challenging 
part is the calculation of the Fock matrix F, which involves 
the core molecular Hamiltonian and the self-consistent po- 
tential. A number of researchers have developed their own 
schemes for performing such ab-initio computations p7|] , 
[ p8[ , p9fl . We accomplish this part by exploiting the fast 
algorithms of Gaussian '98 |3(|, a commercially available 
quantum chemistry software. Aside from the computa- 
tional advantages of using an already well-established soft- 



ware, such a close coupling to the standard tools for ana- 
lyzing molecules makes the chemistry of the system clearer. 
We describe the contacts and the molecule using the sophis- 
ticated LANL2DZ basis set J37J, |Q which incorporates 
relativistic core pseudopotentials. The self-consistent po- 
tential is calculated using DFT with Becke-3 exchange p9| 



and Perdew-Wang 91 correlation |40|. Thus, equipped with 
a SUN workstation, we are able to perform a parameter- 
free analysis of conduction in molecular wires in a few hours 

We illustrate our method using a gold nanowire and a 
phenyl dithiol molecule sandwiched between gold contacts 
and study some previously addressed issues like (1) charge 
transfer and self-consistent band lineup (e.g. jl2)), (2) I-V 
characteristics (e.g. 18 , [ po| , P7| ) and (3) charge den- 
sity and voltage drop (e.g. pq| , p3|). Metallic conduction 
with quantum unit conductance is observed in the gold 
nanowire. Upon reducing the coupling to contacts, the gold 
nanowire exhibits resonant tunneling type of conduction 
just as seen in a phenyl dithiol molecule. The introduction 
of a defect (a stretched bond) in the nanowire gives rise to 
a sharp voltage drop across the impurity as expected. The 
presence of the defect leads to negative-differential resis- 
tance (NDR) in the I-V characteristic of the wire. Periodic 
Friedel oscillations are observed in the charge density near 
the defect, the magnitude of these oscillations decreasing 
as expected upon the introduction of phase-breaking scat- 
tering. 

This paper is organized as follows. Section |n| contains 
a fairly detailed description of the theoretical formulation, 
specifically modeling the influence of the contacts on the 
device subspace through self-energy matrices, and devel- 
oping an appropriate transport formalism to calculate the 
density matrix p (Step 2 above) for the resulting open sys- 
tem under bias. The calculation of the Fock matrix F 
given the density matrix p (Step 1 above) is a standard 
procedure in quantum chemistry and we will not discuss it 
further. Section III shows the results and Section f\/\ briefly 
summarizes the paper. 



II. Theoretical Formulation 
A. Broadening in an open system: Self-energy. 

The concept of self-energy is used in many-body physics 
to describe non-coherent electron-electron and electron- 
phonon interactions. We could do the same in principle 
and use a self-energy function £ p to describe the effect 
of non-coherent interactions of the molecule with its sur- 
roundings (see Fig[l]). For the most part of this paper (one 
exception is made in Fig |ll| a, bottom panel), we will ne- 
glect such non-coherent scattering processes (£ p = 0) be- 
cause the experimental current-voltage characteristics do 
not show any signatures of the molecular vibration spec- 
tra. In addition to £ p , we can use self-energy functions £1 
and £2 to describe the interactions of the molecule with 
the two contacts respectively Q. 

The contact self-energies £1 and £2 arise formally out 
of partitioning an infinite system and projecting out the 



contact Hamiltonians. When an isolated molecule with 
discrete energy levels is contacted to leads to make an infi- 
nite composite system, the energy-dependent one-particle 
retarded Green's function of the complete system is ex- 
pressed in an appropriate basis set as: 



G(E) = [(E + iO+)S - F]~ 



(1) 



where S is the overlap matrix and F is the Fock matrix 
for the whole system. F incorporates the effect of external 
fields, the electronic kinetic energy, electron- nuclear attrac- 
tions, as well as electron-electron interactions, which could 
in principle include Coulomb, exchange and correlation ef- 
fects. The poles of this Green's function lie near the real 
energy axis, and represent the energy levels for the infinite 
system. To extract just the device part of G involving the 
device overlap matrix Sdd and the device Fock matrix Fdd-, 
we utilize the fact that for a matrix 



G = (ES-F)- 1 

ESdd — Fdd 



rt 



D 



G 



d<i 



(2) 



the device part Gdd is given by [ES dd - F dd - tZ>"M) . 
tD^ 1 ^ is the self-energy term that describes the effect of 
the contacts on the device. Only a few surface (interface 
between device and contact) sites give rise to non-zero cou- 
pling elements in r pij , and thus only the surface term of 
.D , the surface Green's function g is needed. We are thus 
left with a reduced device Green's function Gdd given by: 



Gdd(E) = [ESdd - F dd - EiCE) - E 2 (£)]" 



(3) 



where the self-energy matrices Si and £2 are non- 
Hermitian matrices arising from partitioning out contacts 
1 and 2 respectively. Each contact self-energy matrix is re- 
lated to the non-zero part of the corresponding lead-device 
coupling r and the surface Green's function g(E) through: 



El,2 = Tl,2<7l,2Tl,2 



(4) 



The geometry of the bonding between the molecule and 
the contact surface determines the coupling matrices n,2- 
For thiol bonds for example, experimentally it is believed 
that a chemically bonded sulfur atom on a gold surface 
overlaps equally with three gold atoms that form an equi- 
lateral triangle as shown in Fig |l| E|, We use the 
LANL2DZ basis set 0, H to describe both the contacts 
and the molecule. LANL2DZ is a sophisticated basis set 
with relativistic core pseudopotentials and is observed to 
provide a good description of the contact surface density 



A 



B 




Fig. 2 

Central 13 atom part of the 28 atom gold cluster in FCC 

(111) GEOMETRY USED TO CALCULATE THE COUPLING MATRICES S m n 

AND Fmn IN EQ. ^j. A REFERENCE ATOM (ATOM 1) AND ALL IT'S 
NEAREST NEIGHBORS ARE SHOWN IN THIS CENTRAL 13 ATOM PART. 28 
ATOMS ARE USED SO AS TO REDUCE THE 'EDGE' EFFECTS THAT TEND 
TO DESTROY THE SYMMETRIES ASSOCIATED WITH THE FOCK MATRICES 
Fmn (SEE THE DISCUSSION FOLLOWING EQ. [] AND APPENDIX A). 



of states around the Fermi energy 0. The coupling matri- 
ces 71,2 are calculated by using Gaussian '98 to simulate an 
'extended molecule' , Ef| consisting of the molecule un- 
der consideration and an equilateral triangle of three gold 
atoms on either side of the molecule. 

The surface Green's function matrices g\^ are obtained 
recursively for a periodic lattice by the same decimation 
process as for the device. Removing one layer of the contact 
lattice gives back the same surface Green's function, so 
each contact surface Green's function satisfies a recursive 
equation p6[ involving on site and coupling matrices a and 
(3 (diagonal and off-diagonal blocks of ES — F respectively, 
see Fig ||) : 



(5) 



To generalize to a 3-D lead of a specific orientation, we 
follow the procedure in (46) and go to the 2-D fc-space rep- 
resentation for each cross-sectional plane of the lead, so 
that each fc-point effectively acts as an independent 1-D 
problem for which the above recursive formula holds. 

For gold (111) leads, we use Gaussian '98 to extract in- 
plane and out-of-plane nearest-neighbor overlap matrix S 
and Fock matrix F components using a gold cluster (see 



1 For a proper description of the bandstructure of gold away from 
the Fermi energy, we observe that it is necessary to include upto the 
4th or 5th nearest neighboring interactions, because some of the basis 
functions in LANL2DZ are relatively delocalized. 
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Fig ^J) and define Fourier components in the (111) plane: 



71 



ik-(r m —r n ) 



-ik-(r m —r n ) 



(6) 



1-D 
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DEVICE 



where to is an arbitrary gold atom, and n involves a sum 
over to and all its nearest in-plane neighbors, with coor- 
dinate r n . The out-of-plane Fourier components F h ^ and 
S,r are also defined analogously. The F mn and S mn ma- 
trices must all obey the group theoretical symmetry of the 
FCC crystal. This symmetry is satisfied by the S mn matri- 
ces since the overlap between two atoms depends only on 
those two atoms and not on the rest of the atoms in the 
gold cluster. The F mn matrices, however, are obtained via 
a self-consistent calculation (see Fig [TJa) and depend on the 
presence or absence of other atoms in the cluster. Ideally 
we need to simulate an infinite cluster (or crystal) in order 
to get the Fock matrices to obey the symmetry. For practi- 
cal reasons we simulate a finite cluster consisting of 28 gold 
atoms arranged in the FCC (111) geometry (the central 13 
atom part of this cluster is shown in Fig |^) and enforce the 
known symmetry rules for FCC (111) crystal structure on 
the resulting Fock matrices. A detailed description of this 
procedure of enforcing symmetry is given in Appendix A. 
With the correct symmetry imposed, the F j; and S a £ ma- 
trices in Eq. ^| are Hermitian. The gold surface Green's 
function in k space is then obtained by iteratively solving 

''-^-Wl (7) 



where 



and 



(E + M + )S n ^F n 



/% = (£ + *0 + )S bfc - 



The real-space gold (111) surface Green's function matrices 
are then obtained using 



fjn 



1 

N 



J2w 



ik-(r m —r n ) 



(8) 



where N is the number of unit cells in the (111) plane (or 
the number of k points). The surface Green's function ma- 
trices so obtained are independant of the molecule under 
consideration and depend only on the material and geom- 
etry of the contacts. 

We have thus managed to partition the system exactly 
into a device and a lead subspace. The self-energy matrices 
replacing the contacts are non-Hermitian, their real parts 
representing the shift in the molecular energy levels due 
to coupling with the infinite contacts, and their imaginary 
parts representing the broadenings of these levels into a 
continuous density of states. The partitioning makes the 
problem computationally tractable, since the size of the 
self-energy matrices is the same as that of the device Fock 
matrix, even though they represent the effect of infinitely 
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Fig. 3 

Device coupled to a semi-infinite contact (top: 1-D contact, 

BOTTOM: 3-D CONTACT). THE CONTACT SELF-ENERGY DEPENDS ON 
THE DEVICE-CONTACT COUPLING T AND THE CONTACT SURFACE 

Green's function g which is recursively solved for by using 

THE CONTACT ON-SITE AND COUPLING MATRICES a AND f3 (SEE 
TEXT). IN CASE OF A 3-D CONTACT, THE IN-PLANE PERIODICITY IS 
USED TO OBTAIN AN EQUIVALENT 1-D PICTURE WITH K-DEPENDENT 
a k AND j3 k . 



large contacts exactly. In addition, the partitioning de- 
composes the problem into three different subspaces each 
involving a different area of research: (i) the device Fock 
matrix F<id incorporates the quantum chemistry of the in- 
trinsic molecule; (ii) the coupling matrix r involves de- 
tails of the bonding between the molecule and the contact 
(chemisorption, physisorption etc.) and (iii) the surface 
Green's function g involves the surface physics of the metal- 
lic contact, which could in principle be extended to include 
additional effects such as surface states, band-bending, sur- 
face adsorption and surface reconstruction. 

It is desirable to include a few metal atoms as part of 
the device for a number of reasons, (i) The molecule may 
affect a few nearby atoms on the surface of the metallic 
contact. These surface effects will be automatically ac- 
counted for in the self-consistent calculation if a few sur- 
face metal atoms are included as part of the device, (ii) 
Density functional theory is traditionally used for a finite 
system with an integer number of electrons, or periodic 
systems. Extending DFT to a non-neutral open molecular 
subsystem with fractional number of electrons is a topic of 
intense research [|^], [Q. However an extended molecule 
which includes a few metal atoms is effectively charge neu- 
tral and allows the standard DFT formalism to go through, 
(iii) The charges in the molecule are imaged on the metallic 
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Fig. 4 

Comparison of non-self-consistent evaluations of DOS for 
phenyl dithiol with and without a constant self-energy 
approximation. the excellent correspondence is a 
consequence of the relatively flat dos of au(lll) near the 
Fermi energy. 



contact. These image charges need to be taken into account 
in order to accurately calculate the self-consistent poten- 
tial. It is reasonable to expect that the image charges will 
reside on a few surface metal atoms that are close to the 
molecule, and hence the inclusion of these surface atoms 
in the device should account for the image charge effect on 
the self-consistent potential, (iv) Finally, the atoms near 
the two ends of the molecule will have slightly erroneous 
charge densities because the device and contact basis func- 
tions are not orthogonal to each other and the partitioning 
of charge leads to ambiguities at the interfaces (sec Mul- 
liken/Lowdin |49j partitioning). So it is desirable to 'pad' 
the molecule with a few metal atoms on the two ends so as 
to allow an accurate calculation of the molecular charge. 

The results we present in this paper are obtained using 
just the molecule as our 'device' (metal atoms are not in- 
cluded in the device) in order to reduce the computational 
time. Preliminary calculations indicate that including the 
metal atoms improves certain aspects such as the charge 
density on the end atoms. However, such a calculation 
with an extended device takes much more computer time, 
especially with gold contacts. 

The energy dependent self-energy matrices discussed 
above exactly account for the bonding and chemisorption of 
the molecule onto the contact surface, as we discuss in our 
results section. For gold (111) contacts we find that these 
issues are taken care of even if we replace T,(E) with £(£?/) 
(Ef is the gold Fermi energy) as is evident from Fig. ||. This 
helps to reduce the time taken to compute the density ma- 
trix (see Eq. ||). We have developed a fast and elegant 
analytical method to evaluate the density matrix for an 



energy independent self-energy. This method is explained 
in Appendix B. Such a simplification may not be possible 
for platinum contacts having a significant structure in the 
density of states near Ef. This makes the calculation of the 
density matrix computationally quite challenging because 
the correlation function may have sharp peaks in the range 
of integration. Integrating over a complex energy contour 
|42"f , JHofl simplifies the computational complexity to a cer- 
tain degree, but the energy range of integration has to be 
huge in order to take into account all the molecular levels 
(see Appendix B), and a faster scheme is desirable. 

B. Transport: N on- equilibrium Green's function (NEGF) 
formalism. 

The NEGF formalism provides a suitable method for 
calculating the density matrix p for systems under non- 
equilibrium conditions. A tutorial description of this for- 
malism can be found in Here we will simply sum- 
marize the basic relations that can be used (1) to obtain 
p, given the molecular Fock matrix F, the self-energy ma- 
trices Si, £ 2j and the contact electrochemical potentials 
/ii and p,2 and (2) to obtain the electron density n(r) and 
current / from the self-consistent density matrix p. For 
open systems with a continuous density of states, the den- 
sity matrix can be expressed as an energy integral over the 
correlation function —iG < (E) 1 which can be viewed as an 
energy-resolved density matrix: 



p = J dE[-iG < (E)/2Tr} 



The correlation function is obtained from 

-iG< = G(/ 1 r 1 + / 2 r 2 )G+ 



(9) 



(10) 



where fi^(E) are the Fermi functions with electrochemical 
potentials ^i j2 



h,2{E) = 1 + exp 



k B T 



(11) 



G is the Green's function matrix (energy-dependent) in a 
non-orthogonal basis: 



G= (SS-F-Ei-E 2 ) _1 



(12) 



where S is the overlap matrix. The broadening functions 
are the anti-Hermitian components of the self-energy: 



r 1)2 = i[£i, 2 - Si 2 ] 



(13) 



The NEGF equation for p can be interpreted as the fill- 
ing of broadened energy levels (eigenvalues of F) by two 
separate contacts, with GTjC?t representing the density of 
states (DOS) contribution from the ith electrode, and fi 
representing the Fermi function of the ith electrode. Eq. [To| 
assumes that the contacts are reflectionless p4j . 

The converged (see Fig [l] b) density matrix is used to 
obtain the total number as well as the spatial distribution 
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of electrons using 



N 



n{r) 



trace(pS') 



(14) 



where m ,„(r) represent molecular basis functions. The 
converged density matrix may also be used to obtain the 
terminal current H]. For coherent transport f] , we can 
simplify the calculation of the current by using the trans- 
mission formalism where the transmission function ]34| ] : 

T(E) = trace [Tt^G^] 

is used to calculate the terminal current 

/oo 
dET(E) {h{E) - f 2 {E)) 
-OO 

III. Results 



(15) 



A. Equilibrium 

Charge transfer and band lineup. The charge transferred 
between a device and contacts leads to aligning of reference 
energy levels in the two subsystems. Such band lineup is- 
sues are of primary importance in understanding semicon- 
ductor heterojunctions as well as plastic electronics. An 
analogous band-lineup diagram can be obtained on cou- 
pling a molecule to contacts. Fig || shows the surface DOS 
of the contacts and the DOS of the device (a gold nanowire 
in this case). The 6 atom gold chain has discrete energy 
levels that are broadened into a metallic band-like DOS 
on coupling to contacts. For the chain with a broadened 
density of states, one can define a charge-neutrality level 
(CNL) J53|, such that filling up all the states below the 
CNL keeps the device charge-neutral. The Fermi energy 
of gold, approximately -5.1 eV (negative of the bulk gold 
work function) , is noticeably higher than the CNL of a 1-D 
gold nanowire. Owing to this difference electrons flow in 
from the contacts to the device. For a small molecule, the 
capacitative charging energy is large, so the original broad- 
ened levels (dashed) float up self-consistently (solid) due to 
charging till the charge transfer ceases, the band lineup is 
complete and the device and contacts are in chemical equi- 
librium. 

B. N on- equilibrium 

I- V characteristics. A gold nanowire connected to gold 
contacts has a continuous DOS near the equilibrium Fermi 

2 We can modify the transmission formalism to model incoherent 
scattering via Biittiker probes [ pl| . We associate a dephasing strength 
rj and a single electrochemical potential fi p with each Biittiker probe. 
fip is obtained iteratively by asserting that the total current injected 
by the Biittiker probes is zero (f p is the Fermi function with electro- 
chemical potential fi p ): 

(2e/h) J dE Ti p (f p - /i) + J dE T 2p (f p - f 2 ) =0 

This implies that we model incoherent inelastic scattering ^s^opxiosed 
to incoherent elastic scattering which is commonly used |18| , |B3] to 
model phase breaking processes. This procedure was used for the 
bottom panel in Fig. nil a. 
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Fig. 5 

Charge transfer and self-consistent band lineup upon 
coupling a gold nanowire to two gold contacts. the charge 

neutrality level (cnl) of the wire lies below the fermi 
energy of gold, leading to electron flow from contacts to 

the wire. the dos (dashed) of the charge-neutral wire 
redistributes and floats up (solid) self-consistently due to 
charging energies associated with the electron flow. 
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(a) (b) 
Fig. 6 

(a)I-V of a Gold nanowire with six gold atoms forming a 1-D 

CHAIN CONNECTED TO TWO GOLD CONTACTS. An 
ENERGY-INDEPENDENT SELF-ENERGY £(E/) IS USED TO MODEL THE 
CONTACTS. A QUANTUM UNIT CONDUCTANCE OF 2e 2 /h SIGNIFIES 
THAT T,(Ef) ACCOUNTS FOR THE PERFECT TRANSMISSION BY 
SEAMLESSLY JOINING THE NANOWIRE TO THE CONTACTS, (b) I-V OF 
THE SAME GOLD NANOWIRE BUT WITH REDUCED CONTACT COUPLING 
WHICH DECREASES LEVEL BROADENING AND A RESONANT TUNNELING 

TYPE CONDUCTION IS SEEN. THE I-V SHOWS STEP-LIKE BEHAVIOR 
WITH THE CURRENT JUMPING UP WHEN THE CONTACT FERMI LEVELS 
ill AND/OR /J.2 CROSS A MOLECULAR LEVEL (SEE FlGS [j]. ^) . 



\i. = E f - eV/2 



\l 2 = E f + eV/2 



Contact 



Molecule 



Contact 
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Fig. 7 

Mechanism explaining how current flows through a 

MOLECULE COUPLED TO CONTACTS. At EQUILIBRIUM (V = 0) 
flj = i_i 2 = Ef AND NO CURRENT FLOWS. FOR NON-ZERO V THE 

contact Fermi energies separate by an amount eV and a 

SIGNIFICANT CURRENT FLOWS ONLY WHEN A MOLECULAR LEVEL LIES 
IN BETWEEN AND [12 (ALSO SEE FlG ^) . 
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Energy levels of a gold nanowire weakly coupled to 
contacts as a function of applied bias. the i-v for this case 
is shown in Fig |i| b. At zero bias the contact Fermi level is 
CLOSE TO LUMO. AT NON-ZERO BIAS [12 wants to fill up the 
LUMO AND ADD EXTRA ELECTRONS IN THE DEVICE, SO CHARGING 
EFFECTS TEND TO FLOAT THE LEVELS UP TILL fll COMES CLOSE TO 
HOMO AND WANTS TO EMPTY IT. At THIS POINT, MORE ELECTRONS 
ARE LOST THAN GAINED AND CHARGING EFFECTS MAKE THE LEVELS 
FLOAT DOWN PARALLEL TO /ii IN ORDER TO MINIMIZE THE LOSS OF 
ELECTRONS. 



energy (see Fig |^, middle panel), and hence we expect it 
to exhibit metallic conduction as shown in Fig || a. The 
nanowire exhibits a conductance equal to the quantum unit 
conductance of 2e 2 /h w 77 )j,S. This ohmic I-V is a tes- 
timony to the accuracy of our self-energy matrix. Only a 
correct self-energy matrix will get rid of spurious reflections 
at the contact and seamlessly couple the 1-D gold wire with 
the 3-D gold contact. Fig ^ b shows the I-V characteris- 
tic for the same gold nanowire, but with the coupling to 
the contacts reduced by a factor of 4. Due to the reduced 
coupling, the broadening of the levels is much less than the 
separation between the levels, and the wire now exhibits 
resonant conduction |4l]| . There are marked plateaus in 
the I-V curve when a level lies between the contact Fermi 
levels. This is shown schematically in Fig [?]. Fig || shows 
the energy levels of the weakly coupled gold nanowire as a 
function of applied bias. At equilibrium, the contact Fermi 
level lies close to the LUMO level. For a non-zero bias V, 
/ii and /j,2 separate by an amount eV ^ .^2 moves closer to 

3 We choose to split fi± and [12 equally around the equilibrium Fermi 
energy Ef, or fii = Ef — eV/2 and [12 = Ef + eV/2 (see Fig. PI). 
This choice is made in order to be consistent with the Gaussian V8 
convention of applying the linear voltage drop (through the 'field' 
option) symmetrically across the molecule. Specifically, for a molecule 
placed with it's center (along the z-axis) at z = 0, Gaussian '98 



LUMO and wants to fill it up. Charging effects now come 
into the picture and the levels tend to float up and follow 
/12 until /ii comes close to crossing the HOMO level. At 
this point, fi2 wants to put charge into the LUMO while 
/ii wants to empty the HOMO. The number of electrons 
lost is more than the number gained, and now the levels 
tend to go down with [i\ in order to prevent further loss 
of electrons. The current jumps up (see Fig || b) at that 
bias point where fj,2 crosses the LUMO level. Such reso- 
nant conduction is also seen in a phenyl dithiol molecule 
(Fig|l^) with pronounced peaks in the conductance dl/dV 
at resonance. 

The I-V characteristic for a quantum point contact with 
a defect at the center exhibits a weak negative differential 
resistance (NDR) (Fig ||). We artificially stretch a bond 
in the middle of our 6 atom gold chain, which causes the 
left and right local density of states (LDOS) of the chain 
to separately be in equilibrium with the left and right con- 
tacts. Applying a bias then causes the two LDOS to sweep 
past each other within ni. 2 - Owing to the presence of 
some sharp van-Hove singularities in the DOS and some 
smoothened out maxima, there is a progressive alignment 

applies the field such that one end of the molecule is at a voltage 
+V/2 and the other end is at —V/2. One may choose to split )i\ and 
fi2 in some other manner, and the results will not vary provided the 
field is applied consistently. 
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and misalignment of peaks in the LDOS that leads to a 
weak NDR in the I-V. A similar mechanism for NDR has 
been observed in other calculations Q , ||l9[ . 



V(V) 



Fig. 9 

Left and right LDOS for a wire with a defect in the center 

FOR 21 EQUIDISTANT VOLTAGE VALUES FROM TO 2 VOLTS, 
LATERALLY SHIFTED EQUALLY FOR. CLARITY. THE ELECTRONS IN THE 
LEFT AND RIGHT SEGMENTS OF THE WIRE ARE SEPARATELY - IN 
EQUILIBRIUM WITH THE LEFT AND RIGHT CONTACTS, AND THEIR 
RESPECTIVE LDOS FOLLOW THE CORRESPONDING CONTACT 
ELECTROCHEMICAL POTENTIALS (DASHED LINES) BECAUSE OF 
CHARGING EFFECTS SIMILAR TO THOSE SEEN IN FlG ^. As THE TWO 
LDOS SLIDE PAST EACH OTHER WITHIN THE fJ.l-fJ.2 WINDOW (SHOWN 
MAGNIFIED IN THE INSET OF THE BOTTOM FIGURE AT THE FOUR 
VOLTAGE POINTS CIRCLED IN THE I-V GRAPH) THEIR PEAK VALUES 
(TAIL OF THE PEAK FOR THE LEFT LDOS) COME IN AND OUT OF 
RESONANCE, PRODUCING THEREBY' A WEAK NDR IN THE I-V 
CHARACTERISTIC. 





Fig. 10 

i-v of phenyl dithiol connected to two gold contacts. the 
underlying mechanism is resonant tunneling - much like that 
in a gold nanowire weakly coupled to the contacts (see 
Fig [] b). The gap depends on the proximity of the Fermi 
energy to the zero bias homo level if we neglect charging 

effects, while the maximum current at the onset of 
conduction is proportional to the parallel combination of 
broadenings at the molecular energy (see text). for phenyl 
dithiol, the predictions from the above estimates are about 

4 VOLTS AND 40 flA RESPECTIVELY. 



Although most semi-empirical theories qualitatively 
match the shape of the resonant I-V characteristics for 
molecular conductors such as phenyl dithiol (PDT), quan- 
titative agreements between experiment and theory have 
been largely unsatisfactory. The I-V for PDT can be quan- 
tified by the conductance gap and the current value at the 
onset of conduction. If we neglect the charging energy U c 
the gap is roughly given by the proximity of the contact 
Fermi energy to the nearest conducting molecular level, 
Egap ~ 4|.E/ — E mo i \ ||; while the peak current at the on- 
set of conduction is roughly given by 2eT 1 T 2 / (Ti + T 2 ) h, 
evaluated at the energy E mo i . The current onset is smeared 
out over the tail of the molecular level, which in con- 
junction with charging effects leads to a rise in current 
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at the molecular level crossing extended over a voltage 
width AV w (Ti +T2 + Uc) /e. Ab-initio theories are in- 
dispensable for obtaining self-consistently the position of 
the Fermi energy relative to the molecular levels, as well as 
the broadenings of the levels. Both from ab-initio calcula- 
tions for PDT, as well as ab-initio estimates based on the 
above arguments, the predicted conductance gap of PDT, 
set by the molecular HOMO level, comes out to about 4 
Volts, while the maximum current turns out to be tens of 
microamperes. Experimentally both the conductance gap 
and the maximum current are much less than these theo- 
retical estimates. The decrease in current can be attributed 
either to weak coupling between molecular 7r orbitals and 
the s-orbitals of loose gold atoms at the contact surface 
[|54"1| ; WA i or to tunneling between molecular units sepa- 
rately chemisorbed on two ends of a break-junction |p5| . 
Thus a precise experimental knowledge of the contact con- 
ditions is essential prior to appropriate surface modeling. 
The conductance gap is as yet unexplained by the above 
attempts, and requires once again a precise knowledge of 
contact surface conditions which can vary the Fermi energy 
substantially enough to alter the conductance gap. 

Voltage drop and electron density. An applied voltage 
is known to drop largely across the metal-molecule inter- 
face, leading to a weaker drop in the molecule. The precise 
nature of the potential profile is an important input to 
semi-empirical calculations of transport. Tian et al. [ p"8| 
suggested using a flat potential profile inside the molecule, 
with a voltage division factor describing its pos ition. Such 
a flat profile was obtained by Mujica et al. |23|] by solving 
a 1-D Poisson equation, and experimentally measured for 
longer (~ /im) wires by Seshadri and Frisbie How- 
ever, in all these cases the geometry under consideration is 
a series of 2-D charge sheets with potential variations only 
along the wire axis. The 1-D Poisson equation allows vari- 
ations only along one coordinate, while the measurements 
in referred to a self-assembled monolayer (SAM) where 
once again transverse potential variations are screened out 
by the presence of neighboring molecules. In contrast, Lang 
and Avouris p8fl obtained a significant potential drop in a 
carbon atomic wire, which is a consequence of fields pen- 
etrating from transverse directions, as correctly predicted 
for a break-junction geometry by a 3-D Poisson equation. 
Our particular geometry is suited to the break-junction, 
since the Hartree term in Gaussian '98 is calculated for a 
3-D geometry. 

Fig [ll] a shows a gold nanowire weakly coupled to the 
contacts. In Fig[ll] b, the coupling to contacts is strong, but 
a substitutional defect has been incorporated by artificially 
introducing a bond that is longer than the equilibrium gold 
bond length. The voltage drop across the device itself is 
smaller than the applied voltage bias owing to screening 
effects, incorporated self-consistently through the Hartree 
term of our Fock matrix. The plots on top show the net self- 
consistent voltage drop across the wire (solid lines). The 
electrostatic potential does not exhibit peaks at the atomic 
positions because the equilibrium potential profile has been 
subtracted for clarity. Shown for comparison is the linear 




(a) (b) 
Fig. 11 

Electrostatics of (a) a weakly contacted gold wire, and (b) 
a strongly contacted wire with a substitutional defect (a 
stretched bond at the center), in response to a 2 volt 
applied bias across the contacts. top panel: schematic, 
middle panel: voltage drop along the wire (dashed lines 
show a linear drop for comparison), bottom panel: density 

OF ELECTRONS ALONG THE WIRE. In ALL PLOTS, THE CORRESPONDING 

quantities at equilibrium have been subtracted out in order 
to remove peaks near the positions of the nuclei for 
clarity". There is a substantial voltage drop in the wires 
due to transverse fields penetrating the thin (single atom 
cross-section) conductors. the largest potential drop is at 
the barriers (wire-metal interfaces and defects), while the 
electronic charge piles up against the applied bias as 
expected. Superposed on this general polarization of 
charge are friedel oscillations, which die out (dashed line 
in left bottom panel) on increasing incoherent scattering in 
the wire, incorporated through an additional phonon 
self-energy s p introduced through a buttiker probe. 



voltage drop (dashed lines), the solution to Laplace's equa- 
tion in the region. In (b), a large part of the voltage drop 
occurs at the defect as expected. In (a) one might ask why 
the voltage drops across a ballistic device. It is important 
to recognize the distinction between the electrostatic po- 
tential and the electrochemical potential at this point. In 
a ballistic device, it is the electrochemical potential that 
remains constant due to absence of scattering. The elec- 
trostatic potential, however, obeys Poisson's equation and 
may or may not vary depending on the charge density in- 
side the device. For example, a vacuum tube is a ballistic 
device with a linear voltage (electrostatic potential) drop 
due to the absence of any charge. The potential is not well 
screened out due to the thin (one atom) cross-section of the 
wire, which allows transverse fields to penetrate. In con- 
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trast, in a self-assembled monolayer (SAM), the potential 
variation is expected to resemble that predicted in fbsf . 
However, this requires us to disable the Hartree term in 
Gaussian '98 and replace it with our own evaluation (this 
is under current investigation). 

The transferred number of electrons (Fig [H], bottom 
panels) An along the wire obtained from our self-consistent 
solution shows how the voltage drops develop. In both 
plots, there are sizeable oscillations in the charge density, 
which we identify as Friedel oscillations. The equilibrium 
electron density is n — 1/d for the wire (d is the spacing 
between atoms in the wire), leading to a TD Fermi wave- 
vector k,F = mr/2 = ir/2d. The corresponding Friedel os- 
cillation wavelength is 2ir/{2kp) = 2d, about twice an in- 
teratomic spacing, predicting about three oscillation cycles 
in a six atom device, as we indeed observed. Furthermore, 
the amplitudes of the oscillations decrease on incorporat- 
ing an incoherent scattering process in the wire through a 
Biittiker probe (5l) (which shows up as an additional self- 
energy matrix S p representing coupling of the wire with a 
phonon-bath, for example). Superposed on the Friedel os- 
cillations is the general trend of the charge distribution, to 
flow against the applied field (i.e., piling to the left). In (a) 
this piled up electronic charge screens the applied field and 
reduces the slope of the voltage drop in the wire. In (b), 
the charge piles up on each side of the barrier, leading to 
a charge resistivity dipole in the center. This dipole has a 
polar field that is now in the direction of the applied field, 
thus increasing the voltage drop across the barrier . 

IV. Summary 

In this paper we have described a straightforward but 
rigorous procedure for calculating the I-V characteristics 
of molecular wires. The Fock matrix is obtained using 
a Gaussian basis set (LANL2DZ) and the self-consistent 
potential is obtained using density functional theory with 
the B3PW91 approximation. This approach is identical 
to that used in standard quantum chemistry programs 
like Gaussian '98. The difference lies in our use of the 
NEGF formalism to calculate the density matrix from the 
Fock matrix. This allows us to handle open systems far 
from equilibrium which are very different from the isolated 
molecules commonly modeled in quantum chemistry. The 
close coupling with standard quantum chemistry programs 
not only makes the procedure simpler to implement but 
also makes the relation between the I-V characteristics and 
the chemistry of the molecule more obvious. Partitioning 
the contact-molecule-contact system into a molecular sub- 
space and a contact subspace allows us to focus on the 
quantum chemistry of the molecule while exactly taking 
into account the bonding of the molecule with the con- 
tacts, contact surface physics and atomicity etc. We use 
our method to interpolate between two extreme examples 
of transport through a molecular wire connected to gold 
contacts: band conduction in a metallic (gold) nanowire, 
and resonant conduction through broadened, quasidiscrete 
levels of a phenyl dithiol molecule. We examine quanti- 
ties of interest like the I-V characteristic, voltage drop and 



charge density. 
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Appendix A Self-energy: Imposing the FCC 
(III) Symmetry 

In order to calculate the gold (111) surface Green's func- 
tion, we use Gaussian '98 to simulate a 28 atom gold clus- 
ter with FCC (111) geometry and then extract the in-plane 
and out-of-plane overlap and Fock matrix components S mn 
and F mn (see the discussion following Eq. ^ in Section |l|) . 
The cluster, however, does not have the full group theo- 
retical symmetry of the FCC (111) crystal due to edge- 
effects. This symmetry needs to be imposed on the ex- 
tracted Fock matrix elements in order to be consistent with 
the assumptions of two-dimensional periodicity that go into 
our /c-space formalism. The overlap matrix elements auto- 
matically satisfy these symmetry requirements, since they 
involve products of two on-site localized orbitals. 

A typical cluster is shown in Fig with atoms numbered 
from 1 to 13. An FCC (111) surface consists of repetitions 
of such clusters in the sequence ABC ABC ABC . . ., B being 
the plane numbered 1 to 7, and A and C being the planes 
with the rotated triangles above and below it. Each plane is 
parallel to the x-y plane, atom 1 is at the origin and atom 
2 lies on the x-axis. This gives us a z-axis in the (111) 
direction, along which the molecule is assumed to align 
after contacting with the Au surface. The Fock matrix for 
the cluster looks as follows: 



F = 



F\\ F\2 
F21 F22 



F 



1.13 



(16) 



The LANL2DZ basis for gold consists of three sets of 
s, three sets of p and 2 sets of d orbitals, giving us 
3 + 3x3 + 2x5 = 22 orbitals on each gold atom. Each 
of the above matrix components Fn, F12 etc. is therefore 
22 x 22 in size. Using Gaussian '98, we can extract all 
these matrix components separately. For a Au(lll) sur- 
face, however, the elements should be interrelated by the 
point group symmetries of the FCC crystal. Starting with 
two independent matrix components, say, an on-site ele- 
ment Fu that is characteristic of the gold atom, and a 
nearest neighbor element F12 that depends on the lattice 
parameter, we construct the rest of the elements of F using 
the symmetry operations. 

The FCC crystal has the following symmetry: on rotat- 
ing the cluster in Fig. || by 60 degrees about the z-axis 
followed by a reflection in the x-y plane, we get back the 
same cluster. This means there is a unitary matrix M that 
takes F\2 into F±a and so on. The matrix M can be block- 
decomposed into contributions from each of the s, p and d 
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orbital subspaces. Since s orbitals are spherically symmet- 
ric, the corresponding block for the 3 s orbitals is the 3 x 
3 identity matrix. For p orbitals, the M p matrix consists 
of the following reflection (R) and rotation (U) matrices: 





RpUp 


















■(• 


1 












-1 



Un 



cos7r/3 
sin 7r/3 




— sin7r/3 
cos 7r/3 




(17) 



A similar matrix can be constructed for the d-orbital sub- 
space. At the end of the process we thus get a 22 x 
22 M matrix describing rotation-reflection of the entire 
LANL2DZ basis set. 

The nearest neighbor in-plane couplings of atom 1 can 
now be generated by repeated applications of the rotation- 
reflection transformation: 





= M^F 12 M 


Fie 


= M^F 17 M 


Fl 5 


= M^F W M 


Fu 


= M^F 15 M 


F13 


= M'FuM 


F\2 


= M^F 13 M 



(18) 



A simple out-of-plane rotation by 120 degrees can generate 
Fis from F12 ^ . Once again, we explicitly write down 
the contribution for the p orbitals. The d orbitals need 
to be considered analogously while the s orbitals remain 
unchanged. 



Fis — 



M' 

p 




-1/2V3 
5/6 
72/3 




(19) 



All the other matrix components can be similarly obtained, 
using translation and rotation-reflection invariance. For 
example: 



•Fl,12 

-Pi, 10 
F22 
F23 



M^F 1& M 

F n 
F 1A 



(20) 



and so on. 

The crystalline symmetries also constrain the structures 
of the two fundamental matrix components Fn and F\% 

4 For the results shown in this paper, we used the Fig obtained 
from Gaussian '98 instead of generating it from F\2- We thus have 
three independent matrix elements (Fn, F12 and -Fis) instead of±wo 
(-F11 and F±2). Fig has the same structure as that of F12 (see Eq. fejjL . 
Preliminary calculations with _Fis generated from F12 using Eq~|lS| 
do not show a significant change in the results. 



out of which all other components are constructed. -Fn is 
a diagonal block of F and is Hermitian. Fi 2 satisfies the 
following condition: 



(M f ) 3 F12M 3 



F 



(21) 



Three operations of the rotation-reflection transformation 
(equivalently, a coordinate inversion) on F12 generates the 
matrix F15 , which equals F^ by bond- inversion symmetry. 
Translational symmetry implies F51 — F12, which leads to 
the above equation. The form of F12 consistent with the 
above condition is given by: 



F\2 — 




(22) 



where a, e and / are hermitian matrices. The F\2 matrix 
extracted from Gaussian will have this structure for a large 
enough cluster, else we need to impose this condition as 
well. 

With a properly symmetrized set F mn and S mn , the ma- 
trices F a £ and S a ^ in Eq. ^| are Hermitian as expected, and 
we may proceed to calculate the gold (111) surface Green's 
function as outlined in Section Ignoring the symme- 
tries, however, can lead to non-physical results, including 
negative density of states or unsymmetric Green's function 
matrices. 



Appendix B 



Density Matrix: 
Integration 



Analytical 



We have developed a fast and elegant algorithm to eval- 
uate the density matrix given an energy independent self- 
energy. For bulk gold it is known that around the Fermi 
energy the local DOS is approximately a constant (57). It 
is reasonable to expect that all the transport properties of 
the molecule are determined by the molecular levels near 
Ef and the deep levels contribute to the total number of 
electrons by remaining full in the bias range of interest. In 
view of this, We split the energy range of integration in two 
parts and denote the boundary by E},. Ei, is chosen such 
that it is a few electron-volts below the Fermi energy Ef 
and the molecular density of states at Eb is negligible. We 
replace £(-E) by £>(Ef) in the energy range between Ef, to 
Ef. ^| For the energy range from — 00 to Eb we replace 
S(_E) by — it], a constant infinitesimal broadening. The 
only energy dependence in the integrand in Eq. [)] now ap- 
pears through the Green's function G and we can perform 
an analytical integration in an appropriate eigenspace, as 
explained below. 

To evaluate the density matrix, we now need to solve 
integrals of the type (see equations || through [l3| in Sec- 

5 If the DOS is finite at _E{,, there will be a dependence of our results 
on Et, leading to a corresponding ambiguity regarding the precise 
location of Ef at equilibrium. This is an issue that deserves more 
attention, since the precise location of Ef is altered significantly even 
if the total number of electrons, integrated over the entire energy 
range, is miscalculated by a small fraction. As such, the constant self- 
energy approximation used here could also influence Ef in a similar 
way. 
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tion H 



/OO 
dE f(E) G(E) Y G\E) 
-OO 



(23) 



where T is the anti-hermitian part (see Eq. of the en- 
ergy independent self-energy and f(E) is the Fermi func- 
tion with some electrochemical potential /i (see Eq. [ll]) . In 
this paper we assume T — K which modifies the integral 
in the above equation as 



2tt P = dE G(E) Y G\E) 



(24) 



where it is understood that E m i n tends to — oo. Using 
Eq. |l2| for G(E) we may rewrite the above equation as 



2wp = S- 



dE G(E) r g\e) 



S~* 



(25) 



where S is the overlap matrix Y and G(E) are defined as 
T=S-?YS-i 

and 

G(E) ^{EI-Ty 1 
with / being the identity matrix and 

F = 5'-2( j f + s 1 + s 2 )5-3 

In the above equation, Ei 2 are the energy independent 
self-energy matrices as discussed at the beginning of this 
appendix and F is the molecular Fock matrix. 

To solve the integral in Eq. |25| analytically, we work in 
the eigenspace of the energy independent non-Hermitian 
matrix F given by the above equation. Since F is non- 
Hermitian, we need its eigenkets and their dual kets 
{|m)} in order to form a complete basis set. The corre- 
sponding eigenvalues are then {e m } and {e^}. Expanding 
in this mixed basis set, we get 



2np 



= 5-5 



S-i £ dEG(E)(^2\p)(p\ ] j 
YG\e) (X>>v3lJ S ~* 



E„. 



(E - e p ) (E - e* 



(26) 



where we used G{E)\m) = \m)/(E — e m ) and G (E)\rh) = 
\m)/(E-e* m ). 

It is now straightforward to do the above integral an- 
alytically and obtain the density matrix. The current / 
given by Eq. [l5| may also be evaluated analytically using a 
similar procedure. This analytical procedure allows us to 
do the integration extremely fast and avoids errors arising 
out of a numerical integration over a grid. 
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